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Abstract - When a collection of phenotypically diverse organisms compete with each other for 
limited resources, the population can evolve into tightly localised clusters. Past studies have 
neglected the effects of demographic noise and studied the population on a macroscopic scale, 
where cluster formation is found to depend on the shape of the curve describing the decline 
of competition strength with phenotypic distance. Here we show how including the effects of 
demographic noise leads to a radically different conclusion. Two situations are identified: a weak- 
noise regime in which the population exhibits patterns of fluctuation around the macroscopic 
description, and a strong-noise regime where clusters appear spontaneously even in the case that 
all organisms have equal fitness. 
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Introduction. — Competitive interaction between or- 
ganisms has long been recognised as a key component in 
the formation of species and, more widely, entire ecosys- 
tems. A mathematical formulation of this idea was pro- 
vided by MacArthur and Levins [1] over forty years ago 
and variations of their model have been studied ever since, 
with considerable recent interest in the mechanisms by 
which an initially diverse population can aggregate into 
tightly clustered groups [2HE]- In these models, the phe- 
notype of an organism is summarised by a single number 
x, representing a point in the 'niche space' of all possible 
phenotypcs. Writing 3>(x,t) for the size of the popula- 
tion at point x in niche space and t in time, a common 
approach is to study a dynamical system of the 

form: 



d_ 

di 



-i J <f>(x,t)$(y,t)g(x-y)dy. 



(1) 



The three terms on the right of this equation correspond 
to (i) reproduction, where time has been scaled to give this 
rate one everywhere in niche space (ii) mutation, realised 
as diffusion in niche space with strength D, and (iii) death, 
given by the total effect of competition with all other or- 
ganisms. The competition strength between organisms at 



points x and y in niche space is given by the competition 
kernel g(x — y), while the parameter K controls the overall 
carrying capacity of the ecosystem. 

The homogeneous state with total population size K 
is a fixed point of equation ([1]). Starting from this ob- 
servation, a typical analysis considers the effect of small 
perturbations to the homogeneous state [Sj. Often the 
perturbations die out and the population remains evenly 
spread, but in some cases fluctuations of a certain wave- 
length (in niche space) are amplified, eventually leading 
to the formation of sharp peaks separated by empty re- 
gions. Ecologically, these clusters of phenotypically simi- 
lar organisms can be thought of as a prototypical 'species', 
and the model can aid our understanding of the effect of 
competition on species formation. From a physical per- 
spective, the emergence of clusters is a pattern-forming 
transition, and a mathematical analysis reveals that the 
crucial factor is the precise shape of the competition ker- 
nel; species can only form if g has a negative Fourier mode 
and D is sufficiently small. 

One natural choice for the competition kernel is speci- 
fied by the overlap between the organisms' consumption of 
resources, however, it is known that in this case all Fourier 
modes will always be positive, and thus the homogeneous 
state is stable [9]- In general, it may be difficult to find a 
simple ecological argument which justifies the choice of a 
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pattern-forming kernel. Moreover, the common choice of a 
Gaussian kernel is a marginally stable case, meaning that 
a small change to the kernel can lead to very different be- 
haviour. This lack of robustness of the pattern-formation 
mechanism has caused some debate in the literature [6 8 . 

Equation (TT]) is intended to provide a macroscopic sum- 
mary of the collective behaviour of a large number of indi- 
vidual organisms. In reality, the discreteness of the pop- 
ulation means that it will experience some intrinsic de- 
mographic noise; a central assumption of the macroscopic 
theory is that this noise is not important. However, re- 
cent work has revealed that many complex systems ex- 
perience significant noise driven phenomena at scales rel- 
evant to real world observations [T0l[TT]. In this article 
we present the results of a mathematical analysis of noise 
effects in phenotypic competition, with dramatic impli- 
cations for the theoretical understanding of the interplay 
between competition and mutation in ecological systems. 

Our main result is the discovery of two distinct regimes 
of noise- induced phenomena, illustrated in Figure 1. This 
is achieved through the study of an individual-based 
stochastic model of phenotypic competition. We are able 
to show how the macroscopic theory ([1]) may be derived 
directly from the stochastic model (rather than postulated 
a priori), but also how the effects of demographic noise can 
radically alter the observed behaviour. When the strength 
of competition is small in comparison with the rate of mu- 
tation, the model exhibits patterns of fluctuation around 
the homogeneous state. A more extreme effect is found 
when competition strength and mutation rate scale to- 
gether. In this case the population can spontaneously 
cluster into stable macroscopic species, even in the case 
of a completely flat competition kernel, where no selective 
pressure is present at all. 

Stochastic model. — To initiate a proper treatment 
of demographic noise, we must start from a more funda- 
mental description of the system, in terms of the interac- 
tions between individual organisms. We use the following 
microscopic dynamics: 

1. We model niche space as the one-dimensional line seg- 
ment (— 7r,7r]. For mathematical simplicity, and to 
eliminate end effects, we impose periodic boundary 
conditions. 

2. At time t there are N(t) organisms with locations in 
niche space given by x\ , . . . , XjyM . 

3. Each organism reproduces with rate one. The posi- 
tion of the offspring is chosen from a normal distribu- 
tion with variance 2D, centred at the location of the 
parent. 

4. The death rate of organism i is -^J2j9( x i — x j)i 
where, as above, K controls the carrying capacity and 
g is the competition kernel. We take g to be positive, 
symmetric, and normalised so that ^ J g(x)dx — 1. 
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Fig. 1: Noise- induced pattern formation in simulations of phe- 
notypic competition with a unit Gaussian competition kernel 
in (a) the weak-noise regime with mutation rate D = 10~ 3 and 
(b) the strong-noise regime with D = 10 -5 . Dark areas corre- 
spond to a high density of organisms and white space to empty 
regions of niche space. In both cases the macroscopic theory 
predicts that patterns should not form. This, and all other, 
simulations were conducting using the Gillespie algorithm |18] . 
All simulations were initiated with a population of K = 10 3 
organisms distributed uniformly at random throughout niche 
space. To make the images, the locations of the organisms at 
times t = 0, 1, 2, . . . were sorted into 1000 bins, and a grayscale 
value assigned according to the bin occupancy. 

As we will show, equation (fTJ) can be derived directly 
from these dynamics in the limit of large population sizes. 
In that sense, the usual macroscopic theory is a coarse- 
grained description of this microscopic stochastic process. 

A discrete-time version of this model was numerically 
simulated in [12) , and found to exhibit stochastic patterns 
which depend on the shape of the kernel and other pa- 
rameters. We will pursue a mathematical treatment of 
the problem with the advantages that (i) we will be able 
to make concrete, testable, predictions (ii) we can draw 
conclusions about the model on scales which are relevant 
to the real world, and hence too large to be simulated and 
(iii) most importantly, we gain deeper understanding of 
the mechanisms causing the phenomena we observe. 

Expansion of the master equation. At a given 
moment in time, the state of the system is characterised 
by the distribution of individuals 

, N(t) 

= xE<^-^) ■ ( 2 ) 

i=l 

Our mathematical investigation of the model will focus 
on the time evolution of the probability P(<f>, t) of finding 
the system in state <fi at time t. To determine the rate 
of change of P in time, we must consider contributions 
coming from the two processes which alter the system state 
- birth and death. 
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The birth of an organism at location y alters the dis- 
tribution 4> through the addition of another delta func- 
tion. Similarly, the death of organism at y corresponds to 
the removal of the delta peak there. In finite-dimensional 
systems, it is usual to formalise creation and destruction 
processes using so-called step operators (see, for example, 
[15]). We extend that formalism to our case by defining 
the operators A+ and A~ , whose action on a generic func- 
tional F[(j)(x)] is defined to be 



A ± F 

v 



4>{x) 



F 



4>{x) ± j^S(x 



v) 



(3) 



The rate with which new organisms appear at point y is 
the sum of the rates with which each existing organism 
reproduces there. Given that births happen with rate one 
and offspring appear at a Gaussian distance from their 
parent, we determine the birth rate functional to be 

N(t) 

Pt x > <t>) = E r ( x ~ Xi ^> = K I 0(yy( x -y) d y> 

1=1 J 

where we write r{x) for the Gaussian probability density 
function with variance 2D. By a similar logic, the total 
death rate 7 at point x is determined by the product of the 
population density </>(x) with the death rate per organism 
there: 

<P) = 4>{x) ^2 di x ~ x i) = K I <t>( x )<t>(y)9( x -y)dy. 
i=i J 

Combining these two contributions, we obtain the func- 
tional master equation: 



J^P(</>, t) =K J (A- - l) P[ff>, t)p{4>, x) 
+ (A+-l)p(^t) 7 (^x) 



(4) 



To make theoretical progress, we employ a functional vari- 
ant of the Kramers-Moyal expansion [13] , in which the A 
operators are re- written as Taylor series in powers A'" 1 : 

+ <D(K-3), (5) 



A 



1 

1± — 



K 5(f>(x) 2K 2 8<p{xf 



where 8/8<j>{x) denotes functional differentiation. Termi- 
nating the expansion at second order and applying to our 
master equation Q gives the functional Fokker-Planck 
equation 



d_ 

Iff 



P(<f>,t) 



S(j)(x 



1 

2A 



8<f){x) 



■(p(<l>,t)A(<f>,x,y)\ dx dy 

■{p{4>,t)B{4>,x,y)}dxdy, 

(6) 



where 



■A((j), x, y) = (j>(y) r(x - y) - ()>{x)g{x - y) 
B{4>, x , y) = <t>(y) r(x -y) + 4>(x)g(x - y) 



(7) 



This equation will form the basis of our mathematical 
treatment. Aside from its functional form, ([5]) differs from 
a generic Fokker-Planck equation in two regards: firstly, 
the non-local nature of the interactions mean that the drift 
term involves an integral over a second spatial variable; 
secondly, since the step operators in ^ do not appear in 
pairs, the diffusion term is diagonal. 

Depending on the scaling relationship between compe- 
tition strength and mutation rate, the model exhibits two 
different types of behaviour, which we explore in turn. 

Weak-noise: stochastic patterns. — Sending K — > 
00 in ([5]) removes the diffusive part, leaving behind the Li- 
ouvillc equation corresponding to a deterministic variable 
4>{x, t) satisfying 

— cf>(x, t) = tp(y, t) (r(x -y)- 4>(x, t)g{x - y)j dy 

Mr t\ -1- n 

Ox 



d 2 

x ^) + D ^l2^ x ^) ~ I 4>(x,t)cf)(y,t)g(x - y) dy , 



where the second line comes from assuming D is suffi- 
ciently small that only the linear term contributes. Iden- 
tifying $>(x, t) — K<f)(x,t), we obtain precisely the macro- 
scopic model (H}. 

The central limit theorem suggests that the state of the 
microscopic model should be well described by the above 
deterministic equation, plus some stochastic fluctuations 
of order A' -1 / 2 . We look for such fluctuations around the 
homogeneous state 4>{x) = 1/27T by changing variables to 



C(z) 



1 

2^ 



(8) 



whence (O becomes, to leading order in K, 

d ;P(C,t)=- [ [ A(x-y)-^-r{t(y)P(t,t)}dxdy 
8 2 

P((,t)dx, 



dt 



8(( x y 



where A(x) = r(x) — g(x) — 8(x). This equation can be 
solved in Fourier space: writing £„ and A n for the modes 
of £ and A, we find that the £ n are independent (up to 
the condition £_ n = ), zero mean complex Gaussian 
random variables with variances (|Cn| 2 ) evolving according 
to 

l(\ Cn \ 2 ) = 2A n (\( n f) + l. (9) 

Here, and hereafter, angle-brackets refer to the average 
over the noise. Sending t — > 00 here and changing variables 
back to 4>{x) leads to the following description of the long- 
time behaviour of the population density: 



1 1 



n— — oo 



(10) 



To compare our theoretical prediction with the results of 
simulations, it is useful to define a statistic to characterise 
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Fig. 2: Spatial covariance in the weak-noise regime, with pa- 
rameters corresponding to those of Figure 1(a). Blue circles 
show the empirical result approximated by 100 bins and av- 
eraged over 100 independent simulation runs, while the black 
line is the theoretical prediction of equation (|10[) . The peak at 
x = indicates the formation of clusters. 

the structure of the typical states of the system. We choose 
the following scaled measure of spatial covariance: 

H(a;) = 2tt J (tftyMy - x))^ dy , (11) 

where (■ ■ -)oo denotes averaging with respect to the sta- 
tionary (long-time) distribution of <p. This quantity is 
particularly useful as it tells us if the organisms form lo- 
calised clusters or not: a completely flat spatial covariance 
S(x) = 1 would indicate a homogeneous distribution of 
organisms, while a sharp peak at the origin reveals the 
presence of species. In the weak-noise regime, equation 
([IP]) leads to the expression E(a;) = 1 + ± E n (ICn| 2 >e in!E . 
The agreement between simulations and theory is excel- 
lent: sec Figure 2 for an example. 

Many complex systems, particularly those relating to bi- 
ological and chemical processes, exhibit stochastic pattern 
formation in space and / or time pTOlfTT) . This effect can be 
thought of as the stochastic shadow of patterns which oc- 
cur elsewhere in the parameter space of the corresponding 
deterministic system, and it is particularly strong when 
parameters are chosen so that the deterministic system is 
close to the transition between patterned and pattern-free 
phases. This is also true in our case, where the period of 
the stochastic patterns corresponds to the Fourier mode 
which is nearest to becoming unstable. 

As evidenced in Figures 1 and 2, for moderate values of 
K the stochastic patterns can be quite pronounced, mak- 
ing the typical states of the system appear very different to 
the homogeneous distribution predicted by equation (flj. 

Strong-noise: spontaneous species formation. 

Our results in the weak noise regime are based on treating 
the mutation rate D as fixed while taking the limit of large 
carrying capacity K . We will now see that when this as- 
sumption breaks down, the behaviour of the model is rad- 
ically different. In what follows, we will study the special 
case of a flat competition function (g(x) = 1), meaning 
that all organisms experience the same death rate and 
there is no selective pressure at all 1 . Unsurprisingly, the 
macroscopic theory predicts that no species should form in 
this case, that is, the homogeneous distribution is globally 



stable. 

The mathematical study of demographic noise in this 
regime is considerably less straightforward; one can no 
longer appeal to the standard techniques that work for 
weak noise, and must instead develop an approach spe- 
cific to the problem. Our method is to observe that the 
population as a whole, ignoring phenotypic differences, un- 
dergoes logistic growth and should thus rapidly converge 
to a value close to the carrying capacity K. This fact 
may be exploited mathematically through the adiabatic 
elimination of the zeroth Fourier mode of the population 
density. 

As before, we take equation (|6|) for the starting point of 
our analysis, this time expressed in Fourier space. If <j> n 
are the modes of </>, then their joint distribution obeys the 
complex Fokkcr-Planck equation 

n 

(12) 

We are able to isolate (j>o from the rest of the system by 
integrating over the other modes to obtain an equation 
describing noisy logistic growth: 

dP d r 1 

l d 2 (13) 

Ignoring the absorbing state at zero (corresponding to ex- 
tinction), the quasi-stationary distribution of 4>o may be 
found by means of a Wentzcl-Kramcrs-Brillouin expan- 
sion [16]. Briefly, the procedure is as follows: adoption 
of the ansatz P(4>o) = exp{— K S(<fio)} reduces the above 
Fokker-Planck equation to a Hamilton- Jacobi equation for 
S((j>o), plus some correction terms of smaller order in K . 
Solving the Hamilton- Jacobi equation and transforming 
back gives the approximate solution 

P(0o)°c e - 4 ^(i + 27r0o) 4 *. (14) 

For large K, this distribution approaches a peak of width 
approximately (4if7r 2 ) -1 , centred at 1/27T. 

The programme now is to assume that 00 rapidly ap- 
proaches a value close to l/2n, and to approximate the 
behaviour of the other modes based on this premise. The 
most intuitive approach is to work with a stochastic differ- 
ential equation, rather than the Fokkcr-Planck equation. 
Choosing the Ito formalism, equation (fT2"]) (for the flat 
kernel case) is equivalent to the Langevin system 

±6 n = cf> n (e- Dn2 - 2^) + —L^n^t), (15) 



1 We should point out that, as well as being the most important 
case from a paradigmatic point of view, the flat competition kernel 
is a prerequisite for the mathematical techniques we apply. 
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where the r] n (t) are zero mean complex Gaussian white 
noises with covariance structure 

VJtiVmit')) =5(t-t')'/'r l+ rn(e- D( " +m)2 +2^o) • (16) 

The modes must also jointly satisfy the requirement that 
they define the Fourier series of a non-negative function. 
Fortunately, we are able to evaluate the statistics of their 
amplitude without dealing explicitly with this constraint. 
Inserting the approximation <j)o(t) = \j1~K into the n = 
equation yields 

Vo(t) = 0, (17) 
and for the other equations 



dt 



-Dr 



-1 



1 



T)'n(t), (18) 



where the covariance structure of the rj' n (t) is that of the 
original system, conditioned upon the event (|17p . Specifi- 
cally, 



V'n(t)v' m (t)) =5(t-t')i<t> n+m (e 



,-D(r, 



1) 



-Dr, 



1)( 



-Dr, 



A differential equation for the moments (\4> n \ 2 ) can now be 
obtained directly from (|18[) . for example by Ito's formula 
[T7] . Specifically, we find 



dt 



2(e- D 



1 



-Dr, 



1 



(l^l 2 ) 



2A^7T 2 



We are now in a position to examine the behaviour of 
the system under the joint scaling K — > oo and D — > 0, 
where the product of carrying capacity and mutation rate 
approaches a fixed value; we choose KD — > r. In this limit 
we obtain a compact expression for the spatial covariance: 



1 



^ 1 

71 — — OO 



(19) 



This result demonstrates that, in the strong-noise regime, 
the macroscopic equation ([T]) is not a correct description 
of the system as the effects of demographic noise cause 
completely different behaviour. Figure 3 shows an exam- 
ple of the spontaneous emergence of a phenotypic cluster, 
along with comparison between simulations and equation 
([T9"l) for the same parameter values. 

The shape of S(x) for r ranging over several orders of 
magnitude is shown in Figure 4. When r is very small, the 
spatial covariance is sharply peaked at the origin, mean- 
ing that the organisms form a single very tight phenotypic 
cluster. As we move to larger values of r, the typical width 
of this cluster increases; the species is becoming more 
phcnotypically diverse. Eventually as r — > oo we move 
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Fig. 3: Spontaneous species formation with a flat competition 
kernel and KD = 0.01 (K = 10 3 , D = 10~ 5 ). Blue circles show 
the empirical result approximated by 100 bins and averaged 
over 100 independent simulation runs, while the black line is 
the theoretical prediction of equation (|19|l . The strong peak at 
x — indicates the formation of species. 



through the weak-noise regime to a limiting case where 
no species form and the spatial covariance is completely 
flat. Although the calculation presented here applies to 
the case of a flat competition kernel, the effect of sponta- 
neous speciation is robust to changes in the kernel. For 
non-flat kernels it is typical for several distinct species to 
emerge, see Figure 1(b) for an example. 

Discussion. — Simple mathematical models of the 
type inspired by MacArthur and Levins [I] provide a use- 
ful framework in which the effects of phenotypic com- 
petition can be studied. Theoretical work on this topic 
has in the past been limited to the macroscopic level, de- 
scribed by equations such as ([1]). The implicit assumption 
of such work is that the effects of demographic noise are 
not important and may be ignored. In the present analy- 
sis we have taken a ground-up approach to the problem, 
beginning with a description of the system on the funda- 
mental level of interactions between individual organisms 
and culminating in a thorough treatment of the effects of 
demographic noise. In doing so, we have uncovered be- 
haviour which is radically different from that predicted by 
the macroscopic theory. 

Depending on the scaling of mutation rate D and car- 
rying capacity K, two different types of noise effects arc 
present. If D is held constant in the limit K oo the 
typical states of the system can be expressed in terms of 
fluctuations around the homogeneous fixed point of the 
macroscopic equations; we refer to this as the weak-noise 
regime. The precise relationship is given by equation (1101) 
and, as evidenced in Figures 1(a) and 2, the typical states 
of the system are very different to the homogeneous dis- 
tribution predicted by the macroscopic theory. In the 
strong-noise regime, where the mutation rate and carry- 
ing capacity scale together as KD — > r, we find that de- 
mographic noise can become strong enough to cause the 
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Fig. 4: Limiting spatial covariance in the strong-noise regime 
for a range of r values. 



spontaneous formation of well-defined stable species. We 
are able to make theoretical progress by studying this phe- 
nomenon in the paradigmatic case of the flat competition 
kernel, where the limiting form (fTi?)) of the spatial covari- 
ance can be computed. Figure 4 provides an overview of 
the different possible behaviours, interpolating between a 
homogeneous distribution and the formation of tightly lo- 
calised species. Together, these results demonstrate that 
consideration of demographic noise completely overhauls 
our understanding of these models. 

The fact that models of phenotypic competition are ca- 
pable of very different behaviours depending on the rela- 
tionship between mutation rate and carrying capacity may 
have important implications for our understanding of spe- 
ciation in the natural world. Across different taxa, typical 
population sizes vary by many orders of magnitude, whilst 
it is known that there is considerably less variation in the 
rate of genetic mutation [TJ] . Interpreted in the context of 
the present work, our theoretical results suggest that or- 
ganisms belonging to taxa with very large global popula- 
tions (bacteria, foraminifcrs, etc) may form only very loose 
ecological species which could be difficult to tell apart. By 
contrast, the comparatively smaller population sizes com- 
mon to larger creatures imply the spontaneous formation 
of well-defined species. This observation may help explain 
the apparent difficulty in precisely identifying species of 
meiofauna [15], an hypothesis that will receive a careful 
examination in future work. 

The biological interpretation of our results in the con- 
text of species formation is, of course, rather broad since 
the model does not include sexual reproduction, inhomo- 
geneity in niche space or any number of other possible 
complicating factors. Whilst limiting the realism of the 
model, it can be argued that this reductionist approach 
actually strengthens our conclusions; the other factors 
are certainly important in understanding how and where 
species form, but our work shows they cannot be solely 
responsible, since species will form spontaneously even in 
their absence. 

In a wider context, our results highlight the importance 
of theoretical approaches to phenomena of this type. With 



current computing power, individual-based simulations of 
complex systems involving a realistic number of partici- 
pants are not usually possible. It is common to postulate 
that the results of a small simulation can be treated as 
representative of the world at large. However, we have 
shown here that as the number individual agents in a sim- 
ulation increases, radically different behaviours are possi- 
ble depending on how demographic noise interacts with 
the other aspects of a model. This exemplifies the power 
of theoretical methods to provide useful and interesting 
results that could not be obtained from simulations. 
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